require(runjags)
require(coda)
require(lubridate)
require(stringi)

# function to compute a posterior mean and a 95% HPD interval
posterior.summary <- function(x, prob = 0.95) {
  posterior.mean <- mean(x)
  HPD.interval <- HPDinterval(as.mcmc(x), prob = prob)
  c(posterior.mean, HPD.interval)
}

## load the poll data
weekly.PM.data <- read.csv("poll_data.csv")
weekly.PM.data$weekid.2 <- weekly.PM.data$weekid - (weekly.PM.data$pmid - 1) * 233

COVID.data <- read.csv("COVID-19_cases.csv")
COVID.data$date <- ymd(COVID.data$date)
COVID.data$dowid <- wday(COVID.data$date, week_start = 6)
COVID.data$week <- COVID.data$date - COVID.data$dowid + 7

COVID.weekly.mean <- aggregate(COVID.data$case, list(COVID.data$week), mean)
COVID.weekly.mean.diff <- c(rep(0, 200), diff(c(0, COVID.weekly.mean$x))) / 10000
COVID.weekly.mean.ratio <- c(rep(0, 201), diff(COVID.weekly.mean$x) / COVID.weekly.mean$x[-86])

COVID.weekly.Fri <- aggregate(COVID.data$case, list(COVID.data$week),
                              function(x) x[7])
COVID.weekly.Fri$x[1] <- COVID.data$case[6]
COVID.weekly.Fri.diff <- c(rep(0, 200), diff(c(0, COVID.weekly.Fri$x))) / 10000
COVID.weekly.Fri.ratio <- c(rep(0, 203), diff(COVID.weekly.Fri$x)[-1:-2] / COVID.weekly.Fri$x[c(-1:-2, -86)])

# descriptive statistics of the poll data
table(weekly.PM.data$house)
nrow(weekly.PM.data)

## initial values
inits.fun <- function (seed1, seed2, i) {
  set.seed(seed1)
  a1 <- runif(233, 0.4, 0.5)
  a2 <- runif(53, 0.6, 0.7)
  o <- runif(2, 0.01, 0.1)
  b <- runif(10, -0.03, 0.03)
  t <- runif(11, 1, 2)
  list(alpha1 = a1, alpha2 = a2, 
       omega = o, beta.star = b, tau = t, 
       .RNG.name = ifelse(i == 1, 
                          "base::Wichmann-Hill", 
                          ifelse(i == 2, 
                                 "base::Marsaglia-Multicarry", 
                                 "base::Super-Duper")), 
       .RNG.seed = seed2)
}
# preparing three sets of initial values for three MCMC chains
initial.values <- list()
for (i in 1:3) {
  initial.values[[i]] <- inits.fun(i, i, i)
}

#### difference in weekly means: basic ####
result.mean.diff.basic <- run.jags(model = "JAGS_model_1.R", 
                                   monitor = c("gamma", "beta", "tau", "psi", "phi", 
                                               "xi", "omega", "alpha1", "alpha2"), 
                                   data = list(y = weekly.PM.data$pm, 
                                               x = COVID.weekly.mean.diff, 
                                               weekid = weekly.PM.data$weekid.2, 
                                               n.weeks = c(233, 53), 
                                               n.obs = weekly.PM.data$nobs, 
                                               house = weekly.PM.data$house, 
                                               start = c(1, 334), 
                                               end = c(333, 443)), 
                                   inits = initial.values, modules = "glm", 
                                   n.chains = 3, burnin = 1000, sample = 5000, 
                                   adapt = 1000, summarise = FALSE, 
                                   thin = 10, method = "parallel")

# save(result.mean.diff.basic, file = "result_mean_diff_basic.Rdata")
# load("result_mean_diff_basic.Rdata")

gelman.diag(result.mean.diff.basic$mcmc, multivariate = FALSE)
# combining three chains into a single matrix
sims.matrix.mean.diff.basic <- rbind(result.mean.diff.basic$mcmc[[1]], 
                                     result.mean.diff.basic$mcmc[[2]], 
                                     result.mean.diff.basic$mcmc[[3]])
colnames(sims.matrix.mean.diff.basic) <- colnames(result.mean.diff.basic$mcmc[[1]])

# coefficient of the COVID-19 variable
t(round(apply(sims.matrix.mean.diff.basic[, 1:2], 2, 
              function(x) posterior.summary(as.mcmc(x))), 3))

#### figure in the main text ####
## extracting the parameters of cabinet approval
cabinet.approval.samples <- list()
cabinet.approval.samples[[1]] <- 
  sims.matrix.mean.diff.basic[, stri_detect_regex(colnames(sims.matrix.mean.diff.basic), "alpha1")] * 100
cabinet.approval.samples[[2]] <- 
  sims.matrix.mean.diff.basic[, stri_detect_regex(colnames(sims.matrix.mean.diff.basic), "alpha2")] * 100

## computing the posterior means and 95% HPD intervals of cabinet approval ratings
start.date <- ymd(c("2016-4-3", "2020-9-20"))
end.date <- ymd(c("2020-9-13", "2021-8-15"))

cabinet.approval.plot.data <- list()
cabinet.approval.plot.data[[1]] <- data.frame(id = 1:233, 
                                              week = seq(ymd("2016-4-3"), ymd("2020-9-13"), 7), 
                                              rating = colMeans(cabinet.approval.samples[[1]]), 
                                              lower = HPDinterval(as.mcmc(cabinet.approval.samples[[1]]))[, 1], 
                                              upper = HPDinterval(as.mcmc(cabinet.approval.samples[[1]]))[, 2])
cabinet.approval.plot.data[[2]] <- data.frame(id = 234:286, 
                                              week = seq(ymd("2020-9-20"), ymd("2021-9-19"), 7), 
                                              rating = colMeans(cabinet.approval.samples[[2]]), 
                                              lower = HPDinterval(as.mcmc(cabinet.approval.samples[[2]]))[, 1], 
                                              upper = HPDinterval(as.mcmc(cabinet.approval.samples[[2]]))[, 2])

## drawing cabinet approval ratings
png("main_figure.png", width = 90, height = 50, units = "mm", 
    pointsize = 6, res = 1200, type = "cairo")
par(mar = c(3, 3.5, 0, 3.5), lwd = 0.5, family = "Yu Mincho")
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(201, 286), ylim = c(25, 75), 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(cabinet.approval.plot.data[[1]]$id[201:233], 
          cabinet.approval.plot.data[[1]]$id[233:200]), 
        c(cabinet.approval.plot.data[[1]]$lower[201:233], 
          cabinet.approval.plot.data[[1]]$upper[233:200]), 
        border = NA, col = "#B0B0B080")
polygon(c(cabinet.approval.plot.data[[2]]$id, 
          rev(cabinet.approval.plot.data[[2]]$id)), 
        c(cabinet.approval.plot.data[[2]]$lower, 
          rev(cabinet.approval.plot.data[[2]]$upper)), 
        border = NA, col = "#B0B0B080")
points(weekly.PM.data$weekid[weekly.PM.data$weekid > 200], 
       weekly.PM.data$pm[weekly.PM.data$weekid > 200] * 100, 
       pch = 4, cex = 0.4, col = "gray50")
lines(cabinet.approval.plot.data[[1]]$id[201:233], 
      cabinet.approval.plot.data[[1]]$rating[201:233], lwd = 1)
lines(cabinet.approval.plot.data[[2]]$id, 
      cabinet.approval.plot.data[[2]]$rating, lwd = 1)
lines(203:286, 25 + log10(COVID.weekly.mean$x)[-1:-2] * 10)
axis(1, at = c(201, 205, 210, 214, 219, 223, 227, 232, 236, 240, 
               245, 249, 254, 258, 262, 266, 271, 275, 279, 284), 
     labels = NA, lwd = 0.5, family = "serif")
mtext(c("2020", "2021"), side = 1, line = 2, at = c(201, 249), 
      family = "serif")
mtext(c(2:12, 1:9), side = 1, 
      at = c(201, 205, 210, 214, 219, 223, 227, 232, 236, 240, 
             245, 249, 254, 258, 262, 266, 271, 275, 279, 284), 
      line = 1, family = "serif")
mtext(c("��", "�N"), side = 1, line = 1:2, at = 288)
axis(2, at = seq(30, 70, 10), las = 2, lwd = 0.5)
mtext("��\n�t\n�x\n��\n��", side = 2, line = 2.5, las = 2)
axis(4, at = 25 + seq(0, 1, 0.2) * 50, 
     label = c("1", "10", expression(10 ^ 2), expression(10 ^ 3), 
               expression(10 ^ 4), expression(10 ^ 5)), 
     las = 2, lwd = 0.5, family = "serif")
mtext("�T\n��\n��\n�V\n�^\n�R\n��\n�i\n�V\n�K\n��\n��\n��\n��", 
      side = 4, line = 2.5, las = 2)
dev.off()

#### difference in cases on Friday: basic ####
result.Fri.diff.basic <- run.jags(model = "JAGS_model_1.R", 
                                  monitor = c("gamma", "beta", "tau", "psi", "phi", 
                                              "xi", "omega", "alpha1", "alpha2"), 
                                  data = list(y = weekly.PM.data$pm, 
                                              x = COVID.weekly.Fri.diff, 
                                              weekid = weekly.PM.data$weekid.2, 
                                              n.weeks = c(233, 53), 
                                              n.obs = weekly.PM.data$nobs, 
                                              house = weekly.PM.data$house, 
                                              start = c(1, 334), 
                                              end = c(333, 443)), 
                                  inits = initial.values, modules = "glm", 
                                  n.chains = 3, burnin = 1000, sample = 5000, 
                                  adapt = 1000, summarise = FALSE, 
                                  thin = 10, method = "parallel")

# save(result.Fri.diff.basic, file = "result_Fri_diff_basic.Rdata")
# load("result_Fri_diff_basic.Rdata")

gelman.diag(result.Fri.diff.basic$mcmc, multivariate = FALSE)
# combining three chains into a single matrix
sims.matrix.Fri.diff.basic <- rbind(result.Fri.diff.basic$mcmc[[1]], 
                                    result.Fri.diff.basic$mcmc[[2]], 
                                    result.Fri.diff.basic$mcmc[[3]])
colnames(sims.matrix.Fri.diff.basic) <- colnames(result.Fri.diff.basic$mcmc[[1]])

# coefficient of the COVID-19 variable
t(round(apply(sims.matrix.Fri.diff.basic[, 1:2], 2, 
              function(x) posterior.summary(as.mcmc(x))), 3))

#### difference ratio of weekly means: basic ####
result.mean.ratio.basic <- run.jags(model = "JAGS_model_1.R", 
                                    monitor = c("gamma", "beta", "tau", "psi", "phi", 
                                                "xi", "omega", "alpha1", "alpha2"), 
                                    data = list(y = weekly.PM.data$pm, 
                                                x = COVID.weekly.mean.ratio, 
                                                weekid = weekly.PM.data$weekid.2, 
                                                n.weeks = c(233, 53), 
                                                n.obs = weekly.PM.data$nobs, 
                                                house = weekly.PM.data$house, 
                                                start = c(1, 334), 
                                                end = c(333, 443)), 
                                    inits = initial.values, modules = "glm", 
                                    n.chains = 3, burnin = 1000, sample = 5000, 
                                    adapt = 1000, summarise = FALSE, 
                                    thin = 10, method = "parallel")

# save(result.mean.ratio.basic, file = "result_mean_ratio_basic.Rdata")
# load("result_mean_ratio_basic.Rdata")

gelman.diag(result.mean.ratio.basic$mcmc, multivariate = FALSE)
# combining three chains into a single matrix
sims.matrix.mean.ratio.basic <- rbind(result.mean.ratio.basic$mcmc[[1]], 
                                      result.mean.ratio.basic$mcmc[[2]], 
                                      result.mean.ratio.basic$mcmc[[3]])
colnames(sims.matrix.mean.ratio.basic) <- colnames(result.mean.ratio.basic$mcmc[[1]])

# coefficient of the COVID-19 variable
t(round(apply(sims.matrix.mean.ratio.basic[, 1:2], 2, 
              function(x) posterior.summary(as.mcmc(x))), 3))

#### difference ratio of cases on Friday: basic ####
result.Fri.ratio.basic <- run.jags(model = "JAGS_model_1.R", 
                                   monitor = c("gamma", "beta", "tau", "psi", "phi", 
                                               "xi", "omega", "alpha1", "alpha2"), 
                                   data = list(y = weekly.PM.data$pm, 
                                               x = COVID.weekly.Fri.ratio, 
                                               weekid = weekly.PM.data$weekid.2, 
                                               n.weeks = c(233, 53), 
                                               n.obs = weekly.PM.data$nobs, 
                                               house = weekly.PM.data$house, 
                                               start = c(1, 334), 
                                               end = c(333, 443)), 
                                   inits = initial.values, modules = "glm", 
                                   n.chains = 3, burnin = 1000, sample = 5000, 
                                   adapt = 1000, summarise = FALSE, 
                                   thin = 10, method = "parallel")

# save(result.Fri.ratio.basic, file = "result_Fri_ratio_basic.Rdata")
# load("result_Fri_ratio_basic.Rdata")

gelman.diag(result.Fri.ratio.basic$mcmc, multivariate = FALSE)
# combining three chains into a single matrix
sims.matrix.Fri.ratio.basic <- rbind(result.Fri.ratio.basic$mcmc[[1]], 
                                     result.Fri.ratio.basic$mcmc[[2]], 
                                     result.Fri.ratio.basic$mcmc[[3]])
colnames(sims.matrix.Fri.ratio.basic) <- colnames(result.Fri.ratio.basic$mcmc[[1]])

# coefficient of the COVID-19 variable
t(round(apply(sims.matrix.Fri.ratio.basic[, 1:2], 2, 
              function(x) posterior.summary(as.mcmc(x))), 3))

#### difference in weekly means: interaction ####
result.mean.diff.intr <- run.jags(model = "JAGS_model_2.R", 
                                  monitor = c("gamma", "delta", "beta", "tau", "psi", "phi", 
                                              "xi", "omega", "alpha1", "alpha2"), 
                                  data = list(y = weekly.PM.data$pm, 
                                              x = COVID.weekly.mean.diff, 
                                              weekid = weekly.PM.data$weekid.2, 
                                              n.weeks = c(233, 53), 
                                              n.obs = weekly.PM.data$nobs, 
                                              house = weekly.PM.data$house, 
                                              start = c(1, 334), 
                                              end = c(333, 443)), 
                                  inits = initial.values, modules = "glm", 
                                  n.chains = 3, burnin = 1000, sample = 5000, 
                                  adapt = 1000, summarise = FALSE, 
                                  thin = 10, method = "parallel")

# save(result.mean.diff.intr, file = "result_mean_diff_intr.Rdata")
# load("result_mean_diff_intr.Rdata")

gelman.diag(result.mean.diff.intr$mcmc, multivariate = FALSE)
# combining three chains into a single matrix
sims.matrix.mean.diff.intr <- rbind(result.mean.diff.intr$mcmc[[1]], 
                                    result.mean.diff.intr$mcmc[[2]], 
                                    result.mean.diff.intr$mcmc[[3]])
colnames(sims.matrix.mean.diff.intr) <- colnames(result.mean.diff.intr$mcmc[[1]])

# coefficient of the COVID-19 variable * approval rating
t(round(apply(sims.matrix.mean.diff.intr[, 3:4], 2, 
              function(x) posterior.summary(as.mcmc(x))), 3))

#### difference in cases on Friday: interaction ####
result.Fri.diff.intr <- run.jags(model = "JAGS_model_2.R", 
                                 monitor = c("gamma", "delta", "beta", "tau", "psi", "phi", 
                                             "xi", "omega", "alpha1", "alpha2"), 
                                 data = list(y = weekly.PM.data$pm, 
                                             x = COVID.weekly.Fri.diff, 
                                             weekid = weekly.PM.data$weekid.2, 
                                             n.weeks = c(233, 53), 
                                             n.obs = weekly.PM.data$nobs, 
                                             house = weekly.PM.data$house, 
                                             start = c(1, 334), 
                                             end = c(333, 443)), 
                                 inits = initial.values, modules = "glm", 
                                 n.chains = 3, burnin = 1000, sample = 5000, 
                                 adapt = 1000, summarise = FALSE, 
                                 thin = 10, method = "parallel")

# save(result.Fri.diff.intr, file = "result_Fri_diff_intr.Rdata")
# load("result_Fri_diff_intr.Rdata")

gelman.diag(result.Fri.diff.intr$mcmc, multivariate = FALSE)
# combining three chains into a single matrix
sims.matrix.Fri.diff.intr <- rbind(result.Fri.diff.intr$mcmc[[1]], 
                                   result.Fri.diff.intr$mcmc[[2]], 
                                   result.Fri.diff.intr$mcmc[[3]])
colnames(sims.matrix.Fri.diff.intr) <- colnames(result.Fri.diff.intr$mcmc[[1]])

# coefficient of the COVID-19 variable * approval rating
t(round(apply(sims.matrix.Fri.diff.intr[, 3:4], 2, 
              function(x) posterior.summary(as.mcmc(x))), 3))

#### difference ratio of weekly means: interaction ####
result.mean.ratio.intr <- run.jags(model = "JAGS_model_2.R", 
                                   monitor = c("gamma", "delta", "beta", "tau", "psi", "phi", 
                                               "xi", "omega", "alpha1", "alpha2"), 
                                   data = list(y = weekly.PM.data$pm, 
                                               x = COVID.weekly.mean.ratio, 
                                               weekid = weekly.PM.data$weekid.2, 
                                               n.weeks = c(233, 53), 
                                               n.obs = weekly.PM.data$nobs, 
                                               house = weekly.PM.data$house, 
                                               start = c(1, 334), 
                                               end = c(333, 443)), 
                                   inits = initial.values, modules = "glm", 
                                   n.chains = 3, burnin = 1000, sample = 5000, 
                                   adapt = 1000, summarise = FALSE, 
                                   thin = 10, method = "parallel")

# save(result.mean.ratio.intr, file = "result_mean_ratio_intr.Rdata")
# load("result_mean_ratio_intr.Rdata")

gelman.diag(result.mean.ratio.intr$mcmc, multivariate = FALSE)
# combining three chains into a single matrix
sims.matrix.mean.ratio.intr <- rbind(result.mean.ratio.intr$mcmc[[1]], 
                                     result.mean.ratio.intr$mcmc[[2]], 
                                     result.mean.ratio.intr$mcmc[[3]])
colnames(sims.matrix.mean.ratio.intr) <- colnames(result.mean.ratio.intr$mcmc[[1]])

# coefficient of the COVID-19 variable * approval rating
t(round(apply(sims.matrix.mean.ratio.intr[, 3:4], 2, 
              function(x) posterior.summary(as.mcmc(x))), 3))

#### difference ratio of cases on Friday: interaction ####
result.Fri.ratio.intr <- run.jags(model = "JAGS_model_2.R", 
                                  monitor = c("gamma", "delta", "beta", "tau", "psi", "phi", 
                                              "xi", "omega", "alpha1", "alpha2"), 
                                  data = list(y = weekly.PM.data$pm, 
                                              x = COVID.weekly.Fri.ratio, 
                                              weekid = weekly.PM.data$weekid.2, 
                                              n.weeks = c(233, 53), 
                                              n.obs = weekly.PM.data$nobs, 
                                              house = weekly.PM.data$house, 
                                              start = c(1, 334), 
                                              end = c(333, 443)), 
                                  inits = initial.values, modules = "glm", 
                                  n.chains = 3, burnin = 1000, sample = 5000, 
                                  adapt = 1000, summarise = FALSE, 
                                  thin = 10, method = "parallel")

# save(result.Fri.ratio.intr, file = "result_Fri_ratio_intr.Rdata")
# load("result_Fri_ratio_intr.Rdata")

gelman.diag(result.Fri.ratio.intr$mcmc, multivariate = FALSE)
# combining three chains into a single matrix
sims.matrix.Fri.ratio.intr <- rbind(result.Fri.ratio.intr$mcmc[[1]], 
                                    result.Fri.ratio.intr$mcmc[[2]], 
                                    result.Fri.ratio.intr$mcmc[[3]])
colnames(sims.matrix.Fri.ratio.intr) <- colnames(result.Fri.ratio.intr$mcmc[[1]])

# coefficient of the COVID-19 variable * approval rating
t(round(apply(sims.matrix.Fri.ratio.intr[, 3:4], 2, 
              function(x) posterior.summary(as.mcmc(x))), 3))

#### figure in the online appendix ####
# range of each PM's latent approval rating
range(cabinet.approval.plot.data[[1]]$rating)
range(cabinet.approval.plot.data[[2]]$rating)

# estimation results of (gamma + delta * alpha)
mean.diff.Abe.plot <- tcrossprod(sims.matrix.mean.diff.intr[, "gamma[1]"], rep(1, 249)) + 
  tcrossprod(sims.matrix.mean.diff.intr[, "delta[1]"], seq(0.348, 0.596, 0.001))
mean.diff.Suga.plot <- tcrossprod(sims.matrix.mean.diff.intr[, "gamma[2]"], rep(1, 388)) + 
  tcrossprod(sims.matrix.mean.diff.intr[, "delta[2]"], seq(0.293, 0.68, 0.001))
Fri.diff.Abe.plot <- tcrossprod(sims.matrix.Fri.diff.intr[, "gamma[1]"], rep(1, 249)) + 
  tcrossprod(sims.matrix.Fri.diff.intr[, "delta[1]"], seq(0.348, 0.596, 0.001))
Fri.diff.Suga.plot <- tcrossprod(sims.matrix.Fri.diff.intr[, "gamma[2]"], rep(1, 388)) + 
  tcrossprod(sims.matrix.Fri.diff.intr[, "delta[2]"], seq(0.293, 0.68, 0.001))
mean.ratio.Abe.plot <- tcrossprod(sims.matrix.mean.ratio.intr[, "gamma[1]"], rep(1, 249)) + 
  tcrossprod(sims.matrix.mean.ratio.intr[, "delta[1]"], seq(0.348, 0.596, 0.001))
mean.ratio.Suga.plot <- tcrossprod(sims.matrix.mean.ratio.intr[, "gamma[2]"], rep(1, 388)) + 
  tcrossprod(sims.matrix.mean.ratio.intr[, "delta[2]"], seq(0.293, 0.68, 0.001))
Fri.ratio.Abe.plot <- tcrossprod(sims.matrix.Fri.ratio.intr[, "gamma[1]"], rep(1, 249)) + 
  tcrossprod(sims.matrix.Fri.ratio.intr[, "delta[1]"], seq(0.348, 0.596, 0.001))
Fri.ratio.Suga.plot <- tcrossprod(sims.matrix.Fri.ratio.intr[, "gamma[2]"], rep(1, 388)) + 
  tcrossprod(sims.matrix.Fri.ratio.intr[, "delta[2]"], seq(0.293, 0.68, 0.001))

png("appendix_figure.png", width = 5.6, height = 2.8, units = "in", 
    pointsize = 6, res = 1200, type = "cairo")
layout(matrix(1:8, 2, 4, byrow = TRUE))
par(mar = c(3.5, 3.5, 3, 0.5), lwd = 0.5, family = "Yu Mincho")
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.29, 0.68), ylim = c(-4, 2), 
     main = "���{���t�i�����E�T���ρj", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.348, 0.596, 0.001), rev(seq(0.348, 0.596, 0.001))), 
        c(HPDinterval(mcmc(mean.diff.Abe.plot))[, 1], 
          rev(HPDinterval(mcmc(mean.diff.Abe.plot))[, 2])), 
        border = NA, col = "gray80")
lines(seq(0.348, 0.596, 0.001), colMeans(mean.diff.Abe.plot))
abline(h = 0, lty = 3)
axis(1, at = c(0.29, 0.4, 0.5, 0.6, 0.68), lwd = 0.5, family = "serif")
axis(2, lwd = 0.5, family = "serif")
mtext("�O�T�̓��t�x����", side = 1, line = 2.5, cex = 0.8)
mtext("�W\n��", side = 2, line = 2.5, cex = 0.8, las = 2)
par(mar = c(3.5, 3.5, 3, 0.5), lwd = 0.5, family = "Yu Mincho")
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.29, 0.68), ylim = c(-4, 2), 
     main = "���{���t�i�����E���j���j", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.348, 0.596, 0.001), rev(seq(0.348, 0.596, 0.001))), 
        c(HPDinterval(mcmc(Fri.diff.Abe.plot))[, 1], 
          rev(HPDinterval(mcmc(Fri.diff.Abe.plot))[, 2])), 
        border = NA, col = "gray80")
lines(seq(0.348, 0.596, 0.001), colMeans(Fri.diff.Abe.plot))
abline(h = 0, lty = 3)
axis(1, at = c(0.29, 0.4, 0.5, 0.6, 0.68), lwd = 0.5, family = "serif")
axis(2, lwd = 0.5, family = "serif")
mtext("�O�T�̓��t�x����", side = 1, line = 2.5, cex = 0.8)
mtext("�W\n��", side = 2, line = 2.5, cex = 0.8, las = 2)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.29, 0.68), ylim = c(-0.1, 0.06), 
     main = "���{���t�i�ω����E�T���ρj", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.348, 0.596, 0.001), rev(seq(0.348, 0.596, 0.001))), 
        c(HPDinterval(mcmc(mean.ratio.Abe.plot))[, 1], 
          rev(HPDinterval(mcmc(mean.ratio.Abe.plot))[, 2])), 
        border = NA, col = "gray80")
lines(seq(0.348, 0.596, 0.001), colMeans(mean.ratio.Abe.plot))
abline(h = 0, lty = 3)
axis(1, at = c(0.29, 0.4, 0.5, 0.6, 0.68), lwd = 0.5, family = "serif")
axis(2, lwd = 0.5, family = "serif")
mtext("�O�T�̓��t�x����", side = 1, line = 2.5, cex = 0.8)
mtext("�W\n��", side = 2, line = 2.5, cex = 0.8, las = 2)
par(mar = c(3.5, 3.5, 3, 0.5), lwd = 0.5, family = "Yu Mincho")
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.29, 0.68), ylim = c(-0.1, 0.06), 
     main = "���{���t�i�ω����E���j���j", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.348, 0.596, 0.001), rev(seq(0.348, 0.596, 0.001))), 
        c(HPDinterval(mcmc(Fri.ratio.Abe.plot))[, 1], 
          rev(HPDinterval(mcmc(Fri.ratio.Abe.plot))[, 2])), 
        border = NA, col = "gray80")
lines(seq(0.348, 0.596, 0.001), colMeans(Fri.ratio.Abe.plot))
abline(h = 0, lty = 3)
axis(1, at = c(0.29, 0.4, 0.5, 0.6, 0.68), lwd = 0.5, family = "serif")
axis(2, lwd = 0.5, family = "serif")
mtext("�O�T�̓��t�x����", side = 1, line = 2.5, cex = 0.8)
mtext("�W\n��", side = 2, line = 2.5, cex = 0.8, las = 2)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.29, 0.68), ylim = c(-0.4, 0.8), 
     main = "�����t�i�����E�T���ρj", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.293, 0.68, 0.001), rev(seq(0.293, 0.68, 0.001))), 
        c(HPDinterval(mcmc(mean.diff.Suga.plot))[, 1], 
          rev(HPDinterval(mcmc(mean.diff.Suga.plot))[, 2])), 
        border = NA, col = "gray80")
lines(seq(0.293, 0.68, 0.001), colMeans(mean.diff.Suga.plot))
abline(h = 0, lty = 3)
axis(1, at = c(0.29, 0.4, 0.5, 0.6, 0.68), lwd = 0.5, family = "serif")
axis(2, lwd = 0.5, family = "serif")
mtext("�O�T�̓��t�x����", side = 1, line = 2.5, cex = 0.8)
mtext("�W\n��", side = 2, line = 2.5, cex = 0.8, las = 2)
par(mar = c(3.5, 3.5, 3, 0.5), lwd = 0.5, family = "Yu Mincho")
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.29, 0.68), ylim = c(-0.4, 0.8), 
     main = "�����t�i�����E���j���j", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.293, 0.68, 0.001), rev(seq(0.293, 0.68, 0.001))), 
        c(HPDinterval(mcmc(Fri.diff.Suga.plot))[, 1], 
          rev(HPDinterval(mcmc(Fri.diff.Suga.plot))[, 2])), 
        border = NA, col = "gray80")
lines(seq(0.293, 0.68, 0.001), colMeans(Fri.diff.Suga.plot))
abline(h = 0, lty = 3)
axis(1, at = c(0.29, 0.4, 0.5, 0.6, 0.68), lwd = 0.5, family = "serif")
axis(2, lwd = 0.5, family = "serif")
mtext("�O�T�̓��t�x����", side = 1, line = 2.5, cex = 0.8)
mtext("�W\n��", side = 2, line = 2.5, cex = 0.8, las = 2)
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.29, 0.68), ylim = c(-0.06, 0.2), 
     main = "�����t�i�ω����E�T���ρj", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.293, 0.68, 0.001), rev(seq(0.293, 0.68, 0.001))), 
        c(HPDinterval(mcmc(mean.ratio.Suga.plot))[, 1], 
          rev(HPDinterval(mcmc(mean.ratio.Suga.plot))[, 2])), 
        border = NA, col = "gray80")
lines(seq(0.293, 0.68, 0.001), colMeans(mean.ratio.Suga.plot))
abline(h = 0, lty = 3)
axis(1, at = c(0.29, 0.4, 0.5, 0.6, 0.68), lwd = 0.5, family = "serif")
axis(2, lwd = 0.5, family = "serif")
mtext("�O�T�̓��t�x����", side = 1, line = 2.5, cex = 0.8)
mtext("�W\n��", side = 2, line = 2.5, cex = 0.8, las = 2)
par(mar = c(3.5, 3.5, 3, 0.5), lwd = 0.5, family = "Yu Mincho")
plot(NULL, NULL, type = "n", bty = "n", 
     xlim = c(0.29, 0.68), ylim = c(-0.06, 0.2), 
     main = "�����t�i�ω����E���j���j", 
     xlab = "", ylab = "", xaxt = "n", yaxt = "n")
polygon(c(seq(0.293, 0.68, 0.001), rev(seq(0.293, 0.68, 0.001))), 
        c(HPDinterval(mcmc(Fri.ratio.Suga.plot))[, 1], 
          rev(HPDinterval(mcmc(Fri.ratio.Suga.plot))[, 2])), 
        border = NA, col = "gray80")
lines(seq(0.293, 0.68, 0.001), colMeans(Fri.ratio.Suga.plot))
abline(h = 0, lty = 3)
axis(1, at = c(0.29, 0.4, 0.5, 0.6, 0.68), lwd = 0.5, family = "serif")
axis(2, lwd = 0.5, family = "serif")
mtext("�O�T�̓��t�x����", side = 1, line = 2.5, cex = 0.8)
mtext("�W\n��", side = 2, line = 2.5, cex = 0.8, las = 2)
dev.off()